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Using approximate partition functions and a master equation approach, we investigate the statis- 
tical relaxation toward equilibrium in selected CaAr n clusters. The Gaussian theory of absorption 
[J. Chem. Phys., 000, 0000, 2003] is employed to calculate the average photoabsorption intensity 
associated with the 4s 2 — > 4s 1 4p 1 transition of calcium as a function of time during relaxation. In 
CaAr6 and CaAno simple relaxation is observed with a single time scale. CaAri3 exhibits much 
slower dynamics and the relaxation occurs over two distinct time scales. CaAr37 shows much slower 
relaxation with multiple transients, reminiscent of glassy behavior due to competition between dif- 
ferent low-energy structures. We interpret these results in terms of the underlying potential energy 
surfaces for these clusters. 



I. INTRODUCTION 

Clusters of van der Waals atoms or molecules are 
known to exhibit significant finite size effectsEJ in their 
thermodynamic behavior. In particular, their structure 
and low-temperature stable phase can strongly depend 
on size. By doping such clusters with a chromophoric 
atom, spectroscopic techniques are able to characterize 
the structure of the whole complex. Very often, the chro- 
mophore only acts as a small perturbation on the geom- 
etry of the host van der Waals cluster, and spectroscopy 
can then be used as a probe of structure and structural 
changes in the cluster. These ideas, W&ich date back 
to experimental studies in the eighties,0a have received 
some theoretical support in .the recent years Di3 In the 
previous paper in this series,El the influence of tempera- 
ture and cluster size on the photoabsorption spectrum in 
CaAr„ clusters was investigated. In some cases, such as 
CaAri3 or CaAr37, features in the absorption spectrum 
were seen to originate from the competition between spe- 
cific isomers. These results .were obtained using the quan- 
tum superposition method,Q which yields ergodic data by 
construction from a database of local minima on the po- 
tential energy surface (PES). 

However, the equilibrium properties computed from 
the superposition approach obviously provide no infor- 
mation about the way in which equilibrium was estab- 
lished, or how long equilibration would take from a 
given starting distribution. As was shown by Miller 
and coworkers,!!! relaxation can be rather slow when 
the competing minima are separated by large (free) 
energy barriers. In cases such as Ar38, conventional 
molecular dynamics simulations are presently unable to 
reach equilibrium and thus cannot provide estimates 
of the rate constants. However, the master equation 
dynamics approach, introduced in cluster physics by 
Kunz and Berryjiil has been successfully applied to this 
system,L3 although a subsequent discrete path sampling 



study located somewhat faster pathways. Ill In partic- 
ular, it has been used to calculate activation energies 
in the interconversion of (NaCl)J 5 nanocrystalst^l that 
could be compared to those extracted from ion mobility 
measurements!!^ 



It is important to characterize the dynamics of the 
CaAr„ clusters investigated previously at thermodynam- 
ical equilibrium, because the time required for reaching 
equilibrium may not always be accessible in experiments. 
In addition, the most powerful way to investigate phase 
changes through spectroscopy on single clusters may re- 
quire ion trappingO'Ej In principle, this technique can 
be used to examine the system over very long time scales, 
which makes it a possible tool for studying dynamics, in- 
cluding relaxation toward thermal equilibrium, possibly 
via kinetic traps. 

This article is a natural extension to our previous work, 
and provides an investigation of the relaxation dynamics 
in some selected CaAr„ clusters, as evidenced by their 
photoabsorption spectra. In the next section, we briefly 
recall the basic ingredients of the master equation ap- 
proach, and we incorporate quantum corrections to the 
equilibrium probabilities and rate constants. These cor- 
rections are necessary for consistency with the previous 
article.tJ Before actually studying dynamical processes, 
we construct and discuss disconnectivity graphsOc3 in 
Sec. [II. Such graphs are very helpful in elucidating the 
role of the underlying PES on the kinetics. Disconnec- 
tivity graphs are presented and discussed for the clus- 
ters CaArg, CaArio, CaAri3, and CaAr37, which were 
investigated in reference |^. Using the Gaussian theory of 
absorption inspired by Wadi and Pollak,EJt3 we calculate 
the photoabsorption intensity of these clusters as a func- 
tion of both time and excitation energy. The results are 
given and analysed in Sec. |Tv| , before we summarize and 
conclude in Sec. [v| 
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II. MASTER EQUATION DYNAMICS 



In the master equation approach, the interbasin 
dynamics on the energy landscape is described by the 
time evolution of a vector, P(£), whose components are 
the probabilities of the system residing in each of the 
basins at time t. The differential equations governing 
this evolution are 



kjiPiit)} , 



(1) 



where k%j is the rate constant for transitions leading from 
minimum j to minimum i. .A, transition matrix W is 
defined with componentsEZlESt 



where Qj and u>" are the geometric mean normal mode 
frequencies of the minimum and transition state, respec- 
tively. In the low temperature limit (3 — > oo, the rate is 
proportional to T exp[-(Ef -E?)/k B T], where Ef and 
Ej are the energies of the two stationary points includ- 
ing the zero-point energy contributions. While E" — Ej 
is always positive, it may well be that E°"° < E®, leading 
to a divergent rate. In this case the two minima cannot 
be treated independently because of strong quantum de- 
localization effects. However, this situation should not 
arise at high or moderate temperatures where k B T is 
large compared to the zero-point energy, and this is the 
regime of interest in the present work. 
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and Eq. (Q) is solved analytically aftesr-symmetrization 
and diagonalization of the matrix W.li3tL3 The rate con- 
stant kij is the sum over all transition states separating 
minima i and j: 



(3) 



where fc" is given by the usual Rice-Ramsperger-Kassel- 
Marcus theoryE30'E3Eif 00. In the canonical ensemble 
we have 



k?(T) = 
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where k B and h are Boltzmann's and Planck's constants, 
respectively. Z" and Zj are the partition functions corre- 
sponding to the transition state a between i and j, and to 
the minimum j, respectively. Here we use the harmonic 
approximation to model both Z" and Zj , but we incorpo- 
rate quantum corrections by employing the expressions 
for quantum oscillators: 
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In these expressions, hj (resp. h") denotes the order of 
the point group of minimum j (resp. transition state a 
between minima i and j). Ej and E" are the energy of 
these stationary points, and {uijk} and are their 

respective normal mode frequencies, v = 3N — 6 is the 
number of independent degrees of freedom of the cluster, 
and P = l/k B T. In the limit U — > 0, the rate constant 
tends to its usual classical value 
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FIG. 1: Disconnectivity graph of the CaAr6 cluster. The 
energies are in eV. 

Application of the master equation approach requires 
one to build a connected set of minima and the inter- 
well transition statesO To do this we used eigenvector- 
followingpil in conjunction with several methods to sam- 
ple the local minima.E3 For small clusters containing less 
than 11 atoms, systematic searches were started from 
the global minimum, and continued until no new minima 
were found after searching for transition states along ev- 
ery eigenvector of each minimum in the set. For these 
systems the connected set of minima should be nearly 
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complete. For CaAri3, we expect the number of dis- 
tinct local minima to be significantly larger than in 
Arx4, which itself possesses more than 10 4 , excluding 
permutation-inversion isomers. We decided to perform 
only one cycle of systematic eigenvector-following transi- 
tion state searches from each mode of every minimum in 
the database collected for the study at thermodynamic 
equilibrium!] 

The situation becomes quite critical for CaAr 3 7. As 
was shown in our previous studyj3 several minima with 
different structures coexist at low energy, and we were 
not able to build a connected set containing all these 
minima using one systematic eigenvector-following cy- 
cle from the database obtained with parallel temper- 
ing Monte Carlo. For this larger cluster, we selected 
a small set of the low energy minima of interest, 
and we attempted to find discrete pathways (series of 
minima-transition states-minima) connecting them. To 
do this we followed the method described in reference 
|l4| , where a connected path is built up using succes- 
sive double-ended- 



band approacl 



searches. The nudged-elastic- 
was used to generate initial 
guesses for transition states, which were then used as the 
starting points-foc,hybrid eigenvector- following transition 
state searchescjo. These initial paths generally have 
rather high barriers, and so the discrete path sampling 
approach was then employed to provide a more realistic 
account of the rates.li-3 These calculations also produce 
connected databases of minima and transition states that 
should provide a more appropriate picture of the dynam- 
ics, either in visualisations using disconnectivity graphs, 
or from master equation calculations. All the geometry 
optimisations and pathway calcul ati o n s were performed 
with the OPTIM software package JliiKZl for which we im- 
plemented the pairwise potential describing the.-ground 
state potential energy surface of CaAr„ clusters.E3 



III. ENERGY LANDSCAPES 

A convenient way of visualizing complex energy land- 
scape is to construct a disconnectivity graph. Such 
graphs were introduced by Becker and KarplusEj in their 
study of a tetrapeptkLe. They have since been used in a 
much wider context. t3 Disconnectivity graphs show how 
the minima of the PES are connected to each other and, 
more importantly, what is the height of the barriers be- 
tween these minima. Further details about the construc- 
tion of disconnectivity graphs can be found in reference 

H 

In the present work we have included zero-point en- 
ergy corrections for the energies of the minima and tran- 
sition states. Burnham et al. have previously found that 
zero-point effects change the appearance of disconnectiv- 
ity graphs for the water hexamer significantlycil As noted 
in the previous section, the zero-point terms sometimes 
lead to the merging of two (and possibly more) minima. 
However, such events were not observed in CaArg, for 



which we located twelve distinct minima and twenty tran- 
sition states. The disconnectivity graph for this cluster is 
shown in Fig. |l| along with the two lowest energy isomers. 
According to the terminology introduced in reference |3^, 
this graph has a typical 'palm tree' character: the global 
minimum lies at the bottom of a single funnel, and the 
absence of large energy barriers is expected to permit 
relatively fast relaxation toward equilibrium. 
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FIG. 2: Disconnectivity graph of the CaArio cluster. The 
energies are in eV. 

A similar picture was found for the larger cluster, 
CaArio, as shown in Fig. ||. In this case we located 391 
distinct minima connected via 1831 transition states. At 
the bottom of the funnel, the three lowest-energy minima 
lie close together. They are all based on the incomplete 
icosahedra, with different substitional sites for the cal- 
cium atom. As Fig. shows, the two lowest minima are 
extremely close in energy, and the energy barrier that 
separates them is very small (less than 10 -3 eV). Hence 
these two isomers can barely be considered distinct. 

The disconnectivity graph of CaAri3 is shown in Fig. 0. 
The 'palm tree' character of this graph also suggests that 
relatively fast relaxation should be possible to the global 
minimum. However, the bottom of the funnel is slightly 
rugged due to the presence of additional isomers with re- 
spect to what would be expected for Ari4. The four ex- 
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FIG. 3: Disconnectivity graph of the CaAri3 cluster. The right panel shows the lower part of the graph and the secondary 
funnel associated to changing the substitutional site of calcium in the capped icosahedron structure. All energies are in eV. 



tra structures found by changing the substitutional site 
of the calcium atom create a secondary side funnel. They 
are close in energy to the global minimum, and are dis- 
played in Fig. |^ next to a more detailed view of the corre- 
sponding part of the graph. The ordering between these 
isomers is straightforwardly explained by the weaker Ca- 
Ar bond compared to Ar2. The lowest-energy barrier 
between the global minimum and the next four isomers 
lies about 2.19 x 10~ 3 eV/atom above isomer 1, which 
is equivalent to a temperature of about 25 K. Below this 
temperature we expect the relaxation to the global min- 
imum to become significantly slower due to trapping in 
the secondary funnel. 

Finally we consider the larger CaAr37 cluster. Starting 
from the ten lowest minima in the database obtained by 
systematic quenching along a parallel tempering Monte 
Carlo trajectory,tl we were able to establish twenty con- 
nected discrete pathways using the method described in 
reference |l4| New low-energy minima were discovered, 
two of them lower than some in the initial set. Four 
additional pathways were then constructed in the same 
way. As mentioned above, these initial pathways found 
may be far from optimal in terms of rates. The typical 
initial path contained 30-40 minima and the disconnec- 



tivity graph constructed just from the initial pathway 
searches contained 731 distinct minima and 820 transi- 
tion states. It is shown in Fig. [|, together with a selection 
of low-energy structures chosen from the initial sample. 
Clearly this graph does not show the single funnel aspect 
of the previous clusters, and resembles more the 'banyan 
tree' form.E£l 

In fact, the energy landscape of CaAr37 is more com- 
plicated than that of Ar38, which only has two main fun- 
nels associated with truncated octahedral or icosahedral 
structures, respectively. It is closer to the actual multi- 
ple funnel shape p£.the disconnectivity graph in the ionic 
cluster (NaCl) 35-113 This increased complexity is due to 
the larger number of isomers and possible binding sites 
for calcium, and to the existence of a particularly stable 
decahedral isomer. The stabilization of the anti-Mackay 
icosahedral structures (isomer 4j) due to zero-point effects 
is another complicating factor J3 

The CaAr37 cluster exhibits several similarities with 
Ar38, such as the nonicosahedral global minimum and 
the multiple-funnel energy landscape. In thermal equi- 
librium, the low-temperature behavior is essentially gov- 
erned bM the decahedral and Mackay-type icosahedral 
minima.0 The truncated octahedral isomer is not ex- 
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pected to play a major role in the relaxation dynamics, 
except if it contributes significantly to the initial condi- 
tions. 

From Fig. || we see that the decahedral and Mackay 
icosahedra isomers belong to the same superbasin at rel- 
atively low energies. Hence we expect to find signifi- 
cantly faster relaxation than in Ar38, where the trun- 
cated octahedron and the icosahedra are much farther 
apart in configuration space, and separated by a rel- 
atively high barrier. However, to obtain insight into 
the relative stabilities of the isomers at the temperature 
where relaxation is simulated, it may be helpful to-ip-. 
elude entropic effects in the disconnectivity graphs. 
To achieve this goal we simply computed the free energies 
of the minima and transition states instead of potential 
energy minima, replacing Ej (-£■?) by Fj — —ksThiZj 
(Fj* — —ksTlnZj 1 ), where the partition functions were 
obtained in the harmonic approximation, and taken from 
Eqn. §§). 

The free-energy disconnectivity graph of CaAr37 at 
T = 20 K is represented in Fig. |[ As can be clearly 
seen on this figure, the most stable minimum at this tem- 
perature is the anti-Mackay icosahedron isomer (4), and 
the free energy barrier to the decahedral and Mackay- 



icosahedral minima seems rather large. This change in 
the most stable funnel with respect to the T = discon- 
nectivity graph should hinder relaxation, and we there- 
fore anticipate particularly slow dynamics. 

IV. RELAXATION TO EQUILIBRIUM: 
TIME-DEPENDENT ABSORPTION SPECTRA 

We now turn to the relaxation properties of CaAr„ 
clusters resulting from the solution of the master equa- 
tion. Two types of relaxation have been considered, ei- 
ther starting from T = (only the quantum global min- 
imum is occupied) and moving to higher T, or starting 
from a finite temperature equilibrium distribution and 
changing to lower T . The equilibrium distributions were 
calculated within the quantum harmonic superposition 
approach applied to the same database of minima that 
was used to construct the disconnectivity graphs, with- 
out reweighting. For CaAr 3 7, this procedure yields a 
heat capacity peak corresponding to melting located at 
about 25 K, slightly higher than the one obtained in the 
more thorough previous study.B For the other clusters the 
harmonic approximation wthout reweighting performs 
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FIG. 5: Free-energy disconnectivity graph of the CaAr37 clus- 
ter, calculated at T = 20 K using harmonic quantum partition 
functions. The isomers of Fig. W are represented. The energies 
relative to the lowest free-energy minimum are in eV. 

Here we focus on the photoabsorption spectrum as a 
possible probe of cluster structure in a dynamical con- 
text. The absorption intensity of the 4s 2 — > 4s 1 4p 1 ex- 
citation for calcium was calculated using the Gaussian 
theory adapjted and extended from the work of Wadi 
and Pollak.EM The excited state potential energy sur- 
faces were modelled by a diatomics-in-molecules (DIM) 
Hamiltonian, which was fully described in Ref. p8| . 

We first show in Fig. || the time evolution of the ab- 
sorption spectrum of CaArg, starting at t = from the 
zero temperature distribution, and relaxing at T = 25 K, 
slightly above the melting point B The spectrum remains 
essentially doubly-peaked, and at very small times t < 
5 ps the red peak can be separated into two distinct 
lines characteristic of the ground state isomer. The blue 
wing of this peak soon becomes a shoulder. These results 
correspond to fast relaxation, which barely has a spec- 
troscopic signature. Such simple dynamics is perfectly 
consistent with the energy landscape of this cluster, as 
seen from its disconnectivity graph. In the upper panel 
of Fig. ||, the characteristic time constants, r,, of the rate 
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FIG. 6: Photoabsorption spectrum of CaAr6 during the re- 
laxation from the T = distribution toward the T = 25 K 
equilibrium distribution. The discrete spectra [in the quan- 
tum (Q) and classical (C) regimes] of the characteristic times 
of the transition matrix W are represented on a decimal log- 
arithmic scale, on the horizontal axis in the upper part of the 
figure. 



matrix W are represented on a horizontal axis, for both 
the quantum and classical regimes. They are calculated 
from the nonzero eigenvalues, Xi, of W as t; = — 1/Aj. At 
25 K, the largest time constant is a few tens of picosec- 
onds, which is roughly in agreement with the observed 
decay of the spectroscopic signal, although the experi- 
mental relaxationis faster, by about one order of magni- 
tude, than the slowest time constant r max . Such quanti- 
tative discrepancies are typical of the errors that simple 
rate theories based upon harmonic densities of states are 
expected to introduce. Quantum effects are small in this 
system and the relaxation kinetics for CaArio were found 
to be very similar, so detailed discussion is omitted. 

Relaxation from a zero temperature distribution has 
also been investigated in CaAri3 for a temperature jump 
to T = 40 K, which lies close to the melting point of 
this cluster!] At this temperature, the occupation prob- 
ability of the secondary funnel is quite large. The time- 
dependent spectrum in Fig. [7] indeed exhibits a transi- 
tion from the ground state, identified by two absorption 
peaks, to vibrationally excited isomers and the three- 
peak pattern characteristic of calcium located in the 
icosahedral shellJj Monitoring the time-dependent signal 
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FIG. 7: Photoabsorption spectrum of CaAri3 during the re- 
laxation from the T = distribution toward the T = 40 K 
equilibrium distribution. The continuous distributions [in the 
quantum (Q) and classical (C) regimes] of the characteristic 
times of the transition matrix W are represented on a decimal 
logarithmic scale, on the horizontal axis in the upper part of 
the figure. 



shows that relaxation typically occurs on a 50 ps time 
scale, which is again one order of magnitude faster than 
r max in this system, as seen in the upper part of Fig. 0. 

At T = 20 K, the global minimum of CaAri3 is pre- 
dicted to have -an equilibrium occupation probability 
larger than 90%, □ hence the next metastable isomers play 
only a minor role in the thermodynamic properties. Re- 
laxation from a high temperature distribution (50 K) has 
also been simulated using the master equation approach. 
The time-dependent absorption spectrum is plotted in 
Fig. H along with the distribution of characteristic decay 
times from the rate matrix. The variations observed in 
the spectrum are in clear contrast with the previous re- 
sults. During the first 100 ps, the cluster relaxes from 
a broad distribution of isomers towards the main fun- 
nel that corresponds to the capped icosahedron. As in- 
dicated by the two-peak absorption spectrum, the sec- 
ondary funnel has a larger occupation probability about 
5 ns after the dynamics has started. The cluster then fi- 
nally relaxes to its equilibrium distribution following sin- 
gle exponential kinetics, characterized by a rather large 
time constant of about 10 5 ps. This value is still one order 
of magnitude faster than the calculated r max w 10 6 ps, 



represented in the upper part of Fig. |3| The two-step re- 
laxation process observed here can be easily interpreted 
thanks to the disconnectivity graph in Fig. |^. The sec- 
ondary funnel constitutes a kinetic trap, resulting in a 
significant slowing down of the dynamics. While the clus- 
ter easily finds its way into this secondary funnel, escape 
to the global minimum is much slower due to the rather 
large energy barrier. In agreement with previous work, 
this multiple-funnel landscape exhibits a separation of 
timescales in its dynamics. E3o 

We finally discuss the case of CaAr 3 7, for which we 
solved the master equation at T — 20 K, starting from 
the T = initial distribution. The magnitude of the 
characteristic decay times, as computed from the eigen- 
values of the rate matrix, are remarkably large, reach- 
ing values close to one second. The calculated absorp- 
tion spectrum is shown in Fig. |^ as a function of time. 
The choice of a logarithmic scale emphasizes the differ- 
ent time scales present in this system. The multiple fun- 
nels of the energy landscape are in competition below 
20 K, and the cluster undergoes multiple transitions be- 
fore reaching equilibrium. As we deduced from the free 
energy disconnectivity graph in Fig. |^, relaxation from 
to 20 K in this system requires the cluster to escape from 
its initial decahedral/icosahedral funnel. At least four 
significant changes in behaviour are detected through ab- 
sorption spectroscopy, after 10 2 , 10 6 , 10 10 , and 10 12 ps, 
respectively. The transient regime thus appears particu- 
larly slow, in a similar way to the dynamics observed in 
Ar3§.u The spectra in Fig. S show that some intermediate 
minima, which are too high in energy to be significantly 
occupied at equilibrium, may still contribute to the spec- 
trum as transient species with a long lifetime. 

For CaAr37 the set of connected minima was initially 
constructed in a rather artificial manner by means of di- 
rect pathway searches. As mentioned in the previous 
section, these initial paths are usually rather long and 
correspond to high barriers. Therefore the interconver- 
sion rates predicted using this database are probably un- 
realistic and much too low. The discrete path sampling 
tcchniquco was therefore used to optimize the pathway 
between the two lowest decahedral and Mackay icosahe- 
dral minima. We found pathways containing fewer than 
seven minima (including the two end points) with a rate 
constant about eight orders of magnitude larger than that 
of the initial 45-step pathway. This improvement is even 
more impressive than the one found for Ar3sO However, 
in CaAr37 the distance between the two minima in ques- 
tion is rather shorter, because they belong to the same 
funnel. 

We did not attempt to optimize the 23 other initial 
pathways using the discrete path sampling method, since 
the effect of refining the interconversion pathways is not 
expected to affect the complex variation with time of the 
absorption intensity. 
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FIG. 8: Photoabsorption spectrum of CaAri3 during the relaxation from the T = 50 distribution toward the T — 20 K 
equilibrium distribution. The continuous distributions [in the quantum (Q) and classical (C) regimes] of the characteristic 
times of the transition matrix W are represented on a decimal logarithmic scale, on the horizontal axis in the upper part of 
the figure. 



V. CONCLUSIONS 

The complexity of the energy landscapes of CaAr„ 
clusters is reflected not only on their equilibrium prop- 
erties, but also in their relaxation kinetics. By employ- 
ing the master equation approach, quantum corrected 
harmonic partition functions and the Gaussian theory of 
absorption,ll3 we have investigated the time-dependence 
of the photoabsorption intensity. The clusters were cho- 
sen to provide a selection of various finite-size effects. 
Our main result is that spectroscopy may be able to 
probe the isomerization of a given cluster in real time, 
and provide estimates of the interconversion rates them- 
selves. We have found evidence that CaAr„ clusters can 
display simple kinetics (as in CaArg or CaArio), two- 
state kinetics and trapping in auxiliary funnels (aAri3), 
or more intricate kinetics and multiple transient regimes 
(CaAr37), all with a distinct spectroscopic signature, ex- 



periments? 

A number of approximations have been used to achieve 
these results, including the harmonic superposition ap- 
proximation using incomplete samples of minima for the 
larger clusters. The master equation approach also con- 
tains some assumptions, such as local equilibrium in- 
side each basin and Markovian dynamics. Finally, the 
Gaussian theory of absorption also contains some ap- 
proximations. Because of these restrictions, our inves- 
tigation may be semi-quantitative with respect to actual 
CaAr„ clusters, even before we allow for the approximate 
nature of the Hamiltonians used to model -the ground- 
and excited-state potential energy surfaces. E2l Neverthe- 
less, the present study provides the next step beyond the 
analysis of-equilibrium properties discussed in our previ- 
ous paper We believe that it further supports the need 
for future experiments on size-selected, trapped clusters 
studied via spectroscopic techniques. 
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